#ifndef LDPCENT_H_INCLUDED
#define LDPCENT_H_INCLUDED

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Gauss-Markov channel simulator for LDPC codes															      //
// Interlacement of the channel Implementation																	  //
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//#include <windows.h>
//#include <stdio.h>
#include <time.h>
//#include <iostream>
//#include <fstream>
//#include <math.h>
//#include <stdlib.h>
#include <vector>
#include <C:\AMD\Support\LDPC\Headers\archives.h>

#define DEG 4
using namespace std;

int ldpcent(double mi, double pi1, double K, int iErrors, int iRounds, double Eb0, double EbN, double steps[DEG + 1], int iMax, int ene, int ca, double timeT, int degrees[DEG], int threadNumb)
{
//Simulation parameters by user --------------------------------------------------------------------

	unsigned int n=ene,k=ca,itmax=iMax;     // n and k-> matrix sizes   itmax-> maximum number of iterations of Sum-Product Algorithm
	double EbNo=Eb0, EbNoend=EbN, stepsize= steps[0];   //EbNos-> Start and end sizes of noise      stepsize-> size of increments
	double MI=mi, PI1=pi1;							//MI-> memory	PI1-> Probability Bad State
	double Kappa=K;								//Kappa is the impulsive noise

//Simulation parameters by program -----------------------------------------------------------------

	double PI0;										//PI0-> Probability Good State
	double Q, q;									//Transition Probabilities
	unsigned int i,j,indicator=1,iterations=0,ecount1=0,rounds=1,errorcounter1=0;
	unsigned int gaux=0,g=0;
	unsigned int haux=0,h=0;

	double pi=3.14159265358979;
	double Lprod=1,Lsum=0,K1,K2;            //Sum-Product algorithm variables
	double arg1,arg2;
	double U1,U2,V1,V2,W,X,Y,Z;

//Time variables----------------------------------------------------------------

	clock_t sTime, fTime;

//Simulatiom parameters control variables---------------------------------------

	int nonZero = 0;
	int s = 0;
	int degree = 0;

//H loader variables------------------------------------------------------------

	vector<vector<int> >H;                  //Input matrix
	vector<int> aux;

//Codeword, noise and variancy parameters---------------------------------------

	unsigned int state = 0;                 //State of channel
	vector<double> noise(n,0);
	vector<double> r(n,0);
	vector<double> f0(n,0), f1(n,0),A1(n,0),sdVector(n,0);
	double EbNo_Pot,sdg,sdb,sd,BER;
	double R=(double)k/n;

//------------------------------------------------------------------------------

	vector<vector<double> >L1q(n-k,vector<double>(n,0)),L1r(n-k,vector<double>(n,0));
	vector<int> m1(n,0),S1(n-k,0),e1(n,0);   //Message and syndrome vectors

//Reloading last section and grouping simulations-------------------------------

	//Last parameter - 'e'(interlacement), 'g'(genius), '3'(Chase 3 - [genius + memory])
	nonZero = Reload(threadNumb, &iRounds, &iErrors, &Eb0, &timeT, &mi, &pi1, &K, &PI1, &MI, &Kappa, &stepsize, 'e', &s);
	EbNo = Eb0;

//------------------------------------------------------------------------------
//Loading parity check matrix
//------------------------------------------------------------------------------

	H = MatrixLoad(aux, H, n);
	cout<<"Parity check matrices loaded"<<endl;
//------------------------------------------------------------------------------
// Creating Index matrices for storing non zero positions of H matrices
//------------------------------------------------------------------------------
	for(i=0; i<n-k; i++)
	{
		for(j=0; j<n; j++)
		{
			if(H[i][j]==1)
			{
				gaux++;   //finding de 1s in the matrix
			}
			if(g<gaux)
			{
				g=gaux;   //storing the quantity of 1s in g
			}

		}
		gaux=0;
	}

	vector<vector<int> >Index(n-k, vector<int>(g,-1));
	gaux=0;

	for(i=0; i<n-k; i++)
	{
		for(j=0; j<n; j++)
		{
			if(H[i][j]==1)
			{
				Index[i][gaux]=j;       //putting the 1s in position in the Index matrix
				gaux++;                 //in the maximum column index (1s in lines)
			}
		}
		gaux=0;
	}

	for(j=0; j<n; j++)
	{
		for(i=0; i<n-k; i++)
		{
			if(H[i][j]==1)
			{
				haux++;   //finding the 1s in te matrix
			}
			if(h<haux)
			{
				h=haux;   //storing the quantity of 1s em h
			}

		}
		haux=0;
	}

	vector<vector<int> >Indexcol(h,vector<int>(n,-1) );
	haux=0;

	for(j=0; j<n; j++)
	{
		for(i=0; i<n-k; i++)
		{
			if(H[i][j]==1)
			{
				Indexcol[haux][j]=i;    //putting the 1s in position in the Index matrix
				haux++;                 //in the maximum line index (1s in columns)
			}
		}
		haux=0;
	}
	H.clear(); // This is to free the memory used for storing H.

	cout<<"Index matrices loaded"<<endl;
	cout<<"Parity check matrices cleared"<<endl;
	cout<<"Simulation started"<<endl<<endl;

//Simulation------------------------------------------------------------------------------

	EbNo=EbNo-stepsize;

	while(EbNo<EbNoend)
	{
		srand(time(NULL));

// Reinitializing constants---------------------------------------------------

		if (nonZero == 1)			//if is a re-simulation
		{
			rounds=iRounds;
			errorcounter1=iErrors;

            //Making rand() exausted
			int exaust;
			for(int i = 0; i < (4*rounds); i++)	exaust = rand();
		}
		else
		{
			timeT = 0;
			rounds=0;
			errorcounter1=0;
		}
		indicator=1;
		iterations=0;
		ecount1=0;

//Degrees-----------------------------------------------------------------------

		sTime = clock();
		int d = degree;						//actual degree
		if (((rounds >= degrees[0]) || (rounds >= degrees[1]) || (rounds >= degrees[2])) && (nonZero == 0))
		{
			//next degree discovering
			if (rounds < degrees[1])	degree = 1;
			else if (rounds < degrees[2])	degree = 2;
			else if (rounds < degrees[3])	degree = 3;
			else if (rounds >= degrees[3])	degree = 4;

			//walking on steps vector
			if ((d != degree) || (d == 4))	s++;
			if (s <= DEG)	stepsize = steps[s];
		}
		nonZero = 0;
		EbNo=EbNo+stepsize;

//Calculating important values -----------------------------------------------------

		//PI[0] = q/(q+Q);                //PI for good state
		//PI[1] = Q/(q+Q);                //PI for bad state

		PI0 = 1 - PI1;									//PI for good state
		Q = PI1*10000;
		q = PI0*10000;

		EbNo_Pot = pow(10,(EbNo/10));		//Convert EbNo from dB
		sdg = 1.0/(sqrt(2*EbNo_Pot*R));		//Sigma (standard deviation) definition good
		sdb = sqrt(Kappa*sdg*sdg);				//Sigma bad

//Defining initial state -------------------------------------------------------------

		Y = (double)rand()/RAND_MAX;				//Y=[0,1]
		if (Y < PI0)
		{  
			state = 0;								//Define the initial state Good
			sd  = sdg;
		}
		else
		{ 
			state = 1;								//Define the initial state Bad
			sd = sdb;
		}

//Implementing Gauss-Markov channel---------------------------------------------------

		while(errorcounter1<100){ // Stops only when 100 bit errors are found.

			for(i=0;i<n;i++){

//Defining the state--------------------------------------------------------------------

				Z = rand() % (10000 + 1);				//Z=[0,10000]

				if(state == 0){								//for good state
					if(Z <= Q){
						state = 1;							//transition from good to bad
						sd = sdb;
					}
					if(Z > Q){
						state = 0;						//no transition
						sd = sdg;
					}
				}
				else if(state == 1) {										//for bad state
					if(Z <= q){
						state = 0;							//transition from bad to good
						sd = sdg;
					}
					if(Z > q){
						state = 1;						//no transition
						sd = sdb;
					}
				}

				sdVector[i] = sd;					//creating sd vector

//Noise creation ----------------------------------------------------------------

				do
				{
					U1= (double)rand()/RAND_MAX;           // U1=[0,1]
					U2= (double)rand()/RAND_MAX;           // U2=[0,1]
					V1=2*U1 -1;                            // V1=[-1,1]
					V2=2*U2 -1;                            // V2=[-1,1]
					W=V1*V1 + V2*V2;
				} while(W>=1);

				X=sqrt(-2 *log(W) / W) * V1;

				X = sd*X;     //Considering that the noise has zero mean

				noise[i] = X;
			}

			for(i=0;i<n;i++){  // r is the received codeword
				r[i]=1+noise[i];
			}

//Decoding preparation ---------------------------------------------------------


			for(j=0;j<n;j++){ // fi is p[y|x=i]. We consider BPSK modulation, i.e., 0=1 and 1=-1
				f0[j]=(exp((-pow(r[j]-1,2)/(2*sdg*sdg))*PI0)/(sdg*sqrt(2*pi))+exp((-pow(r[j]-1,2)/(2*sdb*sdb))*PI1)/(sdb*sqrt(2*pi)));        //Probability of 0
				f1[j]=(exp((-pow(r[j]+1,2)/(2*sdg*sdg))*PI0)/(sdg*sqrt(2*pi))+exp((-pow(r[j]+1,2)/(2*sdb*sdb))*PI1)/(sdb*sqrt(2*pi)));
			}

			for(i=0;i<n-k;i++){ //Initializing LLR matrix
				for(j=0;j<g;j++){
					if(Index[i][j]!=-1){
						L1q[i][Index[i][j]]=log(f0[Index[i][j]]/f1[Index[i][j]]);
					}
					else{break;}
				}
			}

//Sum Product Algorithm---------------------------------------------------------

			while(indicator==1 && iterations<itmax){

//Calculating the Product-------------------------------------------------------

				for(i=0;i<(n-k);i++){                  //Update rule for check nodes
					for(j=0;j<g;j++){
						if(Index[i][j]!=-1){
							Lprod=Lprod*tanh(L1q[i][Index[i][j]]/2);     //Product
						}
						if(Index[i][j]==-1) break;
					}

					for(j=0;j<g;j++){
						if(Index[i][j]!=-1){ //Note that atan(z)=1/2*[ln(1+z)-ln(1-z)]
							K1= tanh(L1q[i][Index[i][j]]/2);
							K2= (Lprod/K1);
							arg1=1+K2;
							if(arg1<=0){arg1=1.0E-20;}           // This conditional is to avoid underflow
							arg2=1-K2;
							if(arg2<=0){arg2=1.0E-20;}           // This conditional is to avoid underflow
							L1r[i][Index[i][j]]=log(arg1)-log(arg2);     //Filling the matrix to next step
						}
						if(Index[i][j]==-1)break;
					}
					Lprod=1;
				}

//Calculating the Sum-----------------------------------------------------------

				for(j=0;j<n;j++){                 //Update rule for variable nodes
					for(i=0;i<h;i++){
						if(Indexcol[i][j]!=-1){
							Lsum+=L1r[Indexcol[i][j]][j];       //Sum
						}
						if(Indexcol[i][j]==-1)break;
					}

					for(i=0;i<h;i++){
						if(Indexcol[i][j]!=-1){
							L1q[Indexcol[i][j]][j]=log(f0[j]/f1[j])+Lsum-L1r[Indexcol[i][j]][j];
						}
						if(Indexcol[i][j]==-1)break;
					}
					A1[j]=log(f0[j]/f1[j])+Lsum;         //Decision vector
					Lsum=0;
				}

//Filling the message vector---------------------------------------------------

				for(j=0;j<n;j++){              //Termination Rule
					if(A1[j]>=0){m1[j]=0;}
					else m1[j]=1;
				}

//Calculating the Syndrome------------------------------------------------------

				for(i=0;i<n-k;i++){  //Initializing S[]
					S1[i]=0;
				}

				for(i=0;i<n-k;i++){    // Syndrome computation based on the Index matrix S= H*mt
					for(j=0;j<g;j++){
						if(Index[i][j]!=-1){
							S1[i]=(m1[Index[i][j]]+S1[i])%2;
						}
					}
				}


				for(j=0;j<n-k;j++){
					if(S1[j]!=0){
						indicator=1;              //indicates if has Syndrome, if yes, repeats the Sum-Product Algorithm
						break;}
					else indicator=0;         //if don't, step out to errors count
				}

//------------------------------------------------------------------------------

				iterations++;

			}

			indicator=1;
			iterations=0;

			ecount1=0;
			for(i=0; i<n; i++)
			{
				if(m1[i]!=0)
				{
					ecount1++;   //message number of errors count
				}

			}

			if (ecount1>0)   // This indicates error in decoding
			{
				errorcounter1++;
			}

			rounds++;               //number of iterations


			if(rounds%1000==0)
			{
				cout<<"Round number: "<<rounds<<". Detected bit errors: "<<errorcounter1 << " Thread - " << (threadNumb + 1) <<endl;

//Saving Partials ---------------------------------------------------------------

				if(rounds%10000==0)
				{
					fTime = clock();					
					Partials(rounds, errorcounter1, EbNo, fTime, sTime, timeT, MI, PI1, Kappa, stepsize, s, threadNumb);		//Partials function
				}

//-------------------------------------------------------------------------------

			}

		}


		BER=(double)(errorcounter1)/(double)rounds;         //errors number/iterations number

		fTime = clock();

//Saving times -----------------------------------------------------------------

		TimeWrite(EbNo, fTime, sTime, timeT, threadNumb);				//Times print on files

//Printing final results --------------------------------------------------------

		Finals(EbNo, rounds, errorcounter1, BER, threadNumb);			//Finals print on files

		//printing the same on screen
		cout << "Thread - " << (threadNumb + 1) << endl;
		cout<<"EbNo= "<<EbNo<<". Number of rounds for 100 block errors= "<<rounds<<endl;
		cout<<"Total number of errors: "<<errorcounter1<<endl;
		cout<<"BER "<<BER<<endl<<endl;
	}

	return 0;
}

#endif // LDPCENT_H_INCLUDED